class: inverse,left, middle background-image: url(data:image/png;base64,#background.png) background-size: cover <img src="data:image/png;base64,#LOGO_DIPLOMADO.png" width="500px"/> ##Módulo 2: Fundamentos de programación en R para el análisis masivo de datos ### Series de Tiempo Matías Olea <br> <a href="https://orcid.org/0000-0003-0194-7784"> ORCID </a><br> matias.olea@pucv.cl</a><br> .large[<b><a href="https://www.pucv.cl/uuaa/site/edic/base/port/labgrs.html">LabGRS</a> | Octubre 2025</b>] <br> --- class: center,middle background-image: url(data:image/png;base64,#labgrs_logo.png) background-size: 35% --- ## Contenidos .pull-left[ 1) Que es una serie de tiempo. 2) Extracción de series de tiempo (vector) desde raster temporales. 3) Rellenado de datos. 4) Remuestreo temporal. 6) Componentes de una serie y descomposición. 7) Analisis de Anomalias con npphen. ] .pull-right[ <img src="data:image/png;base64,#https://raw.githubusercontent.com/allisonhorst/stats-illustrations/main/rstats-artwork/r_rollercoaster.png" width="650px"/> ] --- ### Series de tiempo Una **serie de tiempo** corresponde a un conjunto de observaciones de variables cuantitativas en el tiempo. Estas variables pueden o no tener distintos comportamiento en el tiempo.
--- ### Series de tiempo Las series las podemos encontrar como estacionarias y no estacionarias .pull-left[ **Estacionaria**. Es aquella cuya media y varianza se mantienen en el tiempo, es decir, son constantes, por lo que facilita más su modelamiento.
] .pull-right[ **No Estacionaria**. Es aquella cuya media y varianza cambian en el tiempo, es decir, no son constantes, por lo que dificulta más su modelamiento.
] --- ### Series de tiempo Las series temporales poseen 4 componentes (3 para algunos autores) los cuales son en simples palabras: -- **Tendencia** (T - trend): corresponde al comportamiento de la media y la varianza en el tiempo. -- **Ciclos** (C - cycle): oscilaciones respecto a la tendencia en largas o medias escalas temporales. Como se desprende de la tendencia, algunos autores no la consideran como un componente propiamente tal. -- **Estacionalidad** (S - seasonality): comportamientos determinados en periodos fijos. Corresponden a cambios dentro de unidades de tiempo de manera regular (como cambios intra-anuales). -- **Aleatorio** (R - random): comportamiento sin un patrón determinado. Es un comportamiento practicamente impredecible. Lo podemos encontrar en la literatura también como "erratico". --- ### Series de tiempo <center><img src="data:image/png;base64,#raster_ts.jpg" width="1000px"/></center> --- ### Series de tiempo #### Lectura de datos ``` r library(terra) library(tidyverse) ``` ``` r setwd("C:/TuDirectorio/") # Leemos fechas desde 2001 y /10000 para pasar de numeros enteros a decimales myd_EVI <- rast("MYD13Q1_EVI_ts491.tif")/10000 # Dividimos por 10000 por escalado del producto tab_dates <- read_csv("MYD13Q1_dates_ts491.csv") # leemos nuestra tabla con fechas ``` ``` r nlyr(myd_EVI) # número de bandas (1 banda = 1 fecha) ``` ``` ## [1] 491 ``` ``` r nrow(tab_dates) # número de fechas ``` ``` ## [1] 491 ``` --- ### Series de tiempo #### Creación serie de tiempo a partir de 1 pixel ``` r xy_plant <- cbind(209482, 6085390) # Coord. Plantación con incendio 2017 xy_nativo <- cbind(209362, 6041365) # Coord. Bosque nativo deciduo ``` ``` r serie_plant <- terra::extract(myd_EVI, xy_plant) %>% as.numeric() # Extraemos serie 1 pix serie_nativo <- terra::extract(myd_EVI, xy_nativo) %>% as.numeric() ``` ``` r tab_pixel <- tibble(dates = rep(tab_dates$date, 2), EVI = c(serie_plant, serie_nativo), Serie = c(rep("Plantacion Mixto", 491), rep("Nativo Deciduo", 491))) ``` ``` r write.csv(tab_pixel, "Tabla_muestreo_pixel.csv", row.names = F) ``` .footnote[ <span style="font-size:9pt"> Como vimos en clases anteriores, en lugar de usar una matriz con las coordenadas XY de un pixel, podemos utilizar un polígono incorporando el estadístico que queremos usar y si queremos incluir los pixeles que toncan los bordes del políno con el argumento touches </span> ] --- ### Series de tiempo #### Creación serie de tiempo a partir de 1 pixel <img src="data:image/png;base64,#GEO455_C8_files/figure-html/unnamed-chunk-13-1.png" width="100%" /> --- ### Series de tiempo #### Rellenado de la serie Existen muchos paquetes que nos permiten realizar rellenos o suavizado de nuestros datos como **gapfill**, **oce**, **imputeTS**, **TSstudio**, o **phenofit**, entre otros. Los más frecuentes para análisis de datos de IV (Índice de Vegetación) proveniente de datos satelitales es **phenopix** o el software **TIMESAT**. Otra opción es la creación de Composiciones (Composites), que realizan agregaciones temporales en ciertos intervalos de tiempo para asegurar una mayor cantidad de datos disponible. Un claro de ejemplo de esto es el producto MOD13Q1 de Modis que realiza composites de 16-días. Para este ejercicio ocuparemos el paquete imputeTS. ``` r library("imputeTS") ``` --- ### Series de tiempo #### Rellenado de datos ``` r tabla_plantacion <- tibble(dates = tab_dates$date, EVI_raw = serie_plant) ``` Para las funciones de imputeTS, necesitamos un objeto clase Time-Series: ``` r ts_plant <- ts(serie_plant, tab_dates$date) ``` ``` r tabla_plantacion$Lin <- na_interpolation(x=ts_plant, # nuestro objeto TS option = "linear", # metodo de rellenado maxgap = Inf) %>% # n de NA consecutivos para rellenar as.numeric() # resultado a vector numerico para incluirlo en la tabla ``` --- ### Series de tiempo #### Remuestreo temporal El remuestreo temporal consiste en realizar **agregaciones** o **agrupaciones** temporales, según una condición de tiempo (mensual, anual, semestral, trimestral, etc.). Comunmente, estos ejercicios se realizan para: - Tener una **periodicidad** o **intérvalos de tiempo regulares** en nuestros datos. Ej: Caso Landsat (16-días 1 satelite / 8-días 2 satelites) - **Disminuir las brechas** por pérdida de datos causados por la presencia de nubes. Ej: Composites de MODIS. De datos diarios a cada 16-dias. - Homologar fuentes de información con **diferentes resoluciones temporales**. Ej: Datos de dendrocronológicos (1-año) con datos satelitales o climáticos. --- class: center,middle background-image: url(data:image/png;base64,#tiempo.png) background-size: 90% background-color: black --- ### Series de tiempo #### Remuestreo temporal Para este ejemplo, utilizaremos nuestra serie de tiempo EVI de MODIS-Aqua de 16-días (23 datos al año) para remuestrearla temporalmene en datos mensuales. **Opción 1: usando tidyverse** ``` r tabla_plantacion$month <- month(tabla_plantacion$dates) # creamos una nueva columna con los meses tabla_plantacion$year <- year(tabla_plantacion$dates) # creamos una nueva columna con los años tabla_mensual <- tabla_plantacion %>% group_by(year, month) %>% # agrupamos por año y mes summarise(EVI_mon = mean(Lin, na.rm = T)) # agrupamos usando la media tabla_mensual$dates <- paste(tabla_mensual$year, tabla_mensual$month,"15", sep = "-") %>% as.Date() # construimos vector de fechas para las agrupaciones ``` --- ### Series de tiempo #### Remuestreo temporal <img src="data:image/png;base64,#GEO455_C8_files/figure-html/unnamed-chunk-19-1.png" width="100%" /> --- ### Series de tiempo #### Remuestreo temporal **Opción 2: usando xts** ``` r library(xts) # la versión para raster es "rts" xts_plantacion <- xts(tabla_plantacion$Lin, # vector numerico con la variable tabla_plantacion$dates) # vector de fechas EVI_mensual <- apply.monthly(xts_plantacion, FUN = mean) # agregación mensual ``` --- ### Series de tiempo #### Remuestreo temporal <img src="data:image/png;base64,#GEO455_C8_files/figure-html/unnamed-chunk-21-1.png" width="100%" /> --- ### Series de tiempo #### Test de estacionariedad ``` r library("tseries") ``` **Kwiatkowski-Phillips-Schmidt-Shin (KPSS) test** **H<sub>0</sub>**: La serie de tiempo tiene una tendencia estacionaria. **H<sub>A</sub>**: La serie de tiempo tiene una tendencia no estacionaria. --- ### Series de tiempo #### Test de estacionariedad ``` r kpss.test(tabla_plantacion$Lin, null = "Trend", lshort = F) # Plantacion Mixto (rellenada) ``` ``` ## ## KPSS Test for Trend Stationarity ## ## data: tabla_plantacion$Lin ## KPSS Trend = 0.30546, Truncation lag parameter = 17, p-value = 0.01 ``` Se rechaza H<sub>0</sub> por lo tanto la serie es no estacionaria. <img src="data:image/png;base64,#GEO455_C8_files/figure-html/unnamed-chunk-24-1.png" width="100%" /> --- ### Series de tiempo #### Test de estacionariedad ``` r kpss.test(na.omit(serie_nativo), null = "Trend", lshort = F) # Nativo (con NA) ``` ``` ## ## KPSS Test for Trend Stationarity ## ## data: na.omit(serie_nativo) ## KPSS Trend = 0.10429, Truncation lag parameter = 17, p-value = 0.1 ``` No se rechaza H<sub>0</sub> por lo tanto la serie es estacionaria. <img src="data:image/png;base64,#GEO455_C8_files/figure-html/unnamed-chunk-26-1.png" width="100%" /> --- ### Series de tiempo #### Descomposición de la serie ``` r ts_mensual <- ts(tabla_mensual$EVI_mon, # vector con la variable start = c(2002, 7), # año y mes de inicio de la serie end = c(2023, 10), # año y mes de termino de la serie frequency = 12) # frecuencia de datos al año dec_plant <- decompose(ts_mensual) # descomposición ``` --- ### Series de tiempo <img src="data:image/png;base64,#GEO455_C8_files/figure-html/unnamed-chunk-29-1.png" width="100%" /> --- ### Series de tiempo: Analisis con npphen <center><img src="data:image/png;base64,#Screenshot_2.png" height="500px" /></center> --- ### Series de tiempo: Analisis con npphen <center><img src="data:image/png;base64,#Screenshot_1.png" height="500px" /></center> --- ### Series de tiempo: Analisis con npphen <center><img src="data:image/png;base64,#Screenshot_3.png" height="500px" /></center> --- ### Series de tiempo: Analisis con npphen <center><img src="data:image/png;base64,#Screenshot_4.png" height="500px" /></center> --- ### Series de tiempo: Analisis con npphen <center><img src="data:image/png;base64,#Screenshot_5.png" height="500px" /></center> --- ### Series de tiempo: Analisis con npphen <center><img src="data:image/png;base64,#Screenshot_6.png" height="500px" /></center> --- ### Series de tiempo: Analisis con npphen <center><img src="data:image/png;base64,#Screenshot_7.png" height="500px" /></center> --- ### Series de tiempo: Analisis con npphen <center><img src="data:image/png;base64,#Screenshot_8.png" height="500px" /></center> --- ### Series de tiempo: Analisis con npphen <center><img src="data:image/png;base64,#Screenshot_9.png" height="500px" /></center> --- ### Series de tiempo: Analisis con npphen <center><img src="data:image/png;base64,#Screenshot_10.png" height="500px" /></center> --- ### Series de tiempo: Analisis con npphen <center><img src="data:image/png;base64,#Screenshot_11.png" height="500px" /></center> --- ### Series de tiempo: Analisis con npphen <center><img src="data:image/png;base64,#Screenshot_12.png" height="500px" /></center> --- ### Series de tiempo: Analisis con npphen <center><img src="data:image/png;base64,#Screenshot_13.png" height="500px" /></center> --- ### Series de tiempo: Construcción de fenología Vamos a definir el periodo de referencia o estable como las fechas entre el 2010 y 2016 ``` r inicio <- which(tabla_plantacion$year == 2010) %>% first() termino <- which(tabla_plantacion$year == 2016) %>% last() evi_referencia <- tabla_plantacion$EVI_raw[inicio:termino] date_referencia <- tabla_plantacion$dates[inicio:termino] ``` --- ### Series de tiempo: Construcción de fenología ``` r library(npphen) PhenKplot(x = evi_referencia, dates = date_referencia, h = 2, # periodo h = 2 Julio - Junio; h = 1 = Enero - Diciembre xlab ='Dia del ciclo', ylab ='EVI', rge = c(0, 0.6)) # rge = rango teorico ``` <img src="data:image/png;base64,#GEO455_C8_files/figure-html/unnamed-chunk-31-1.png" width="100%" /> --- ### Series de tiempo: Calculo de anomalias ``` r anom <- ExtremeAnom(x = tabla_plantacion$EVI_raw*10000, # en valores ENTEROS dates = tabla_plantacion$dates, h =2, refp = c(inicio:termino), # periodo de referencia anop = c(1:491), # periodo al que se calcula anomalias rge = c(0, 6000), # rango teorico en ENTEROS output = 'anomalies') ``` --- ### Series de tiempo: Calculo de anomalias ``` r barplot(anom / 10000, # para re-escalar los valores enteros col = ifelse(anom > 0, 'green', 'red'), names = '', main = 'Anomalias') ``` <img src="data:image/png;base64,#GEO455_C8_files/figure-html/unnamed-chunk-33-1.png" width="100%" /> --- ### Series de tiempo: Calculo de anomalias extremas ``` r anom_extrema <- ExtremeAnom(x = tabla_plantacion$EVI_raw*10000, # en valores ENTEROS dates = tabla_plantacion$dates, h =2, refp = c(inicio:termino), # periodo de referencia anop = c(1:491), # periodo al que se calcula anomalias rge = c(0, 6000), # rango teorico en ENTEROS output = 'clean', # filtradas por solo extremas rfd = 0.95) # frecuencia definida extrema ``` --- ### Series de tiempo: Calculo de anomalias ``` r barplot(anom_extrema / 10000, # para re-escalar los valores enteros col = ifelse(anom_extrema > 0, 'green', 'red'), names = '', main = 'Anomalias') ``` <img src="data:image/png;base64,#GEO455_C8_files/figure-html/unnamed-chunk-35-1.png" width="100%" /> --- ### Bibliografía (2012) E. B. Brooks, V. A. Thomas, R. H. Wynne and J. W. Coulston, "Fitting the Multitemporal Curve: A Fourier Series Approach to the Missing Data Problem in Remote Sensing Analysis," in IEEE Transactions on Geoscience and Remote Sensing, vol. 50, no. 9, pp. 3340-3353. doi: 10.1109/TGRS.2012.2183137. (2004) Jin Chen, Per. Jönsson, Masayuki Tamura, Zhihui Gu, Bunkei Matsushita, Lars Eklundh, A simple method for reconstructing a high-quality NDVI time-series data set based on the Savitzky–Golay filter, Remote Sensing of Environment, Volume 91, Issues 3–4, (2022) Kong, D., McVicar, T. R., Xiao, M., Zhang, Y., Peña-Arancibia, J. L., Filippa, G., Xie, Y., Gu, X. phenofit: An R package for extracting vegetation phenology from time series remote sensing. Methods in Ecology and Evolution, 13, 1508– 1527. https://doi-org.pucv.idm.oclc.org/10.1111/2041-210X.13870 (2015) G. Yang, H. Shen, L. Zhang, Z. He and X. Li, "A Moving Weighted Harmonic Analysis Method for Reconstructing High-Quality SPOT VEGETATION NDVI Time-Series Data," in IEEE Transactions on Geoscience and Remote Sensing, vol. 53, no. 11, pp. 6008-6021. doi: 10.1109/TGRS.2015.2431315. --- class: inverse middle 